/********************************************************************************
 *    Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH    *
 *                                                                              *
 *              This software is distributed under the terms of the             * 
 *              GNU Lesser General Public Licence (LGPL) version 3,             *  
 *                  copied verbatim in the file "LICENSE"                       *
 ********************************************************************************/
#include "EtfNeutronWall.h"
#include "EtfNeutronWallFirstHits.h"
#include "EtfNeutronWallPoint.h"
#include "EtfNeutronWallGeo.h"
#include "EtfNeutronWallGeoPar.h"

#include "FairVolume.h"
#include "FairGeoVolume.h"
#include "FairGeoNode.h"
#include "FairRootManager.h"
#include "FairGeoLoader.h"
#include "FairGeoInterface.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"
#include "EtfDetectorList.h"
#include "EtfStack.h"

#include "TClonesArray.h"
#include "TVirtualMC.h"
#include "TGeoManager.h"
#include "TGeoBBox.h"
#include "TGeoCompositeShape.h"
#include "TGeoTube.h"
#include "TGeoMaterial.h"
#include "TGeoMedium.h"

#include <iostream>
using std::cout;
using std::endl;

EtfNeutronWall::EtfNeutronWall()
  : FairDetector("EtfNeutronWall", kTRUE, kEtfNeutronWall)
{
    ResetParameters();
    fEtfNeutronWallPointCollection=new TClonesArray("EtfNeutronWallPoint");
    flGeoPar = new TList();
    flGeoPar->SetName( GetName());
    fEtfNeutronWallFirstHits = new TClonesArray("EtfNeutronWallFirstHits");
}

EtfNeutronWall::EtfNeutronWall(const char* name, Bool_t active)
  : FairDetector(name, active, kEtfNeutronWall)
{
    ResetParameters();
    fEtfNeutronWallPointCollection=new TClonesArray("EtfNeutronWallPoint");
    flGeoPar = new TList();
    flGeoPar->SetName( GetName());
    fEtfNeutronWallFirstHits = new TClonesArray("EtfNeutronWallFirstHits");
}



EtfNeutronWall::~EtfNeutronWall()
{
  if ( flGeoPar ) delete flGeoPar;
  if (fEtfNeutronWallPointCollection) {
    fEtfNeutronWallPointCollection->Delete();
    delete fEtfNeutronWallPointCollection;
  
  if (fEtfNeutronWallFirstHits) {
    fEtfNeutronWallFirstHits->Delete();
    delete fEtfNeutronWallFirstHits;
  }
  
  }
}

void EtfNeutronWall::Initialize()
{
  /**
   * WORKAROUND needed for Geant4 in MT mode
   * Call AddSensitiveVolume for sensitive volumes in order to fill
   * thread-local FairModule::svList.
   */
  DefineSensitiveVolumes();
  FairDetector::Initialize();
  FairRuntimeDb* rtdb= FairRun::Instance()->GetRuntimeDb();
  EtfNeutronWallGeoPar* par=(EtfNeutronWallGeoPar*)(rtdb->getContainer("EtfNeutronWallGeoPar"));
    // Initialise variables from Birk law
  Double_t dP = 1.032 ;
  // Set constants for Birk's Law implentation
  fBirkC0 =  1.;
  fBirkC1 =  0.013/dP;
  fBirkC2 =  9.6e-6/(dP * dP);
}

Bool_t  EtfNeutronWall::ProcessHits(FairVolume* vol)
{
  /** This method is called from the MC stepping */
  Getpripoint();
  Int_t cp1=-1;
  Int_t cp2=-1;
  Int_t volId1=-1;
  Int_t volId2=-1;
  volId1 =  gMC->CurrentVolID(cp1);
  volId2 =  gMC->CurrentVolOffID(1, cp2);
  //Set parameters at entrance of volume. Reset ELoss.
  if ( TVirtualMC::GetMC()->IsTrackEntering() ) {
    fELoss  = 0.;
    fTime   = TVirtualMC::GetMC()->TrackTime() * 1.0e09;
    fLength = TVirtualMC::GetMC()->TrackLength();
    TVirtualMC::GetMC()->TrackPosition(fPosIn);
    TVirtualMC::GetMC()->TrackMomentum(fMomIn);
  }

  // Sum energy loss for all steps in the active volume
  fELoss += TVirtualMC::GetMC()->Edep();
  Double_t lightYield =  gMC->Edep() ;
  
    // Correction for all charge states
  if (gMC->TrackCharge()!=0) {
    Double_t birkC1Mod = 0;
     // Apply correction for higher charge states
    if (fBirkC0==1){
      if (TMath::Abs(gMC->TrackCharge())>=2)
        birkC1Mod=fBirkC1*7.2/12.6;
      else
        birkC1Mod=fBirkC1;
    }

    Double_t dedxcm=0.;
    if (gMC->TrackStep()>0)
    {
      dedxcm=1000.*gMC->Edep()/gMC->TrackStep();
//      cout<<"Edep: "<<gMC->Edep()<<" Trackstep: "<<gMC->TrackStep()<<endl;
      lightYield=lightYield/(1.+birkC1Mod*dedxcm+fBirkC2*dedxcm*dedxcm);
//      cout<<"dedxcm "<<dedxcm<<" factor "<<1./(1.+birkC1Mod*dedxcm+fBirkC2*dedxcm*dedxcm)
//      <<" light yield "<<lightYield<<endl;
      fLightYield=fLightYield+lightYield;
//      cout<<"fLightYield: "<<fLightYield<<endl;
    }
  }

  // Create EtfNeutronWallPoint at exit of active volume
  if ( TVirtualMC::GetMC()->IsTrackExiting()    ||
       TVirtualMC::GetMC()->IsTrackStop()       ||
       TVirtualMC::GetMC()->IsTrackDisappeared()   ) {
    fTrackID  = TVirtualMC::GetMC()->GetStack()->GetCurrentTrackNumber();
    fVolumeID = vol->getMCid();
    if (fELoss == 0. ) { return kFALSE; }
    gMC->TrackPosition(fPosOut);
    gMC->TrackMomentum(fMomOut);
    
    if (gMC->IsTrackExiting()) {
      const Double_t* oldpos;
      const Double_t* olddirection;
      Double_t newpos[3];
      Double_t newdirection[3];
      Double_t safety;
    
      gGeoManager->FindNode(fPosOut.X(),fPosOut.Y(),fPosOut.Z());
      oldpos = gGeoManager->GetCurrentPoint();
      olddirection = gGeoManager->GetCurrentDirection();
    
  //     cout << "1st direction: " << olddirection[0] << "," << olddirection[1] << "," << olddirection[2] << endl;

      for (Int_t i=0; i<3; i++){
        newdirection[i] = -1*olddirection[i];
      }
    
      gGeoManager->SetCurrentDirection(newdirection);
     //  TGeoNode *bla = gGeoManager->FindNextBoundary(2);
      safety = gGeoManager->GetSafeDistance();


      gGeoManager->SetCurrentDirection(-newdirection[0],-newdirection[1],-newdirection[2]);
    
      for (Int_t i=0; i<3; i++){
        newpos[i] = oldpos[i] - (3*safety*olddirection[i]);
      }

      fPosOut.SetX(newpos[0]);
      fPosOut.SetY(newpos[1]);
      fPosOut.SetZ(newpos[2]);
    }
    
    AddHit(fTrackID, fVolumeID, cp2,cp1,TVector3(fPosIn.X(),  fPosIn.Y(),  fPosIn.Z()),
           TVector3(fMomIn.Px(), fMomIn.Py(), fMomIn.Pz()),TVector3(fPosOut.X(),  fPosOut.Y(),  fPosOut.Z()),
           TVector3(fMomOut.Px(), fMomOut.Py(), fMomOut.Pz()), fTime, fLength,
           fELoss);

    // Increment number of EtfNeutronWall det points in TParticle
    EtfStack* stack = (EtfStack*) TVirtualMC::GetMC()->GetStack();
    stack->AddPoint(kEtfNeutronWall);
  }

  return kTRUE;
}

void EtfNeutronWall::Register()
{

  /** This will create a branch in the output tree called
      EtfNeutronWallPoint, setting the last parameter to kFALSE means:
      this collection will not be written to the file, it will exist
      only during the simulation.
  */

  if ( ! gMC->IsMT() ) {
    FairRootManager::Instance()->Register("EtfNeutronWallPoint", "EtfNeutronWall",
                                          fEtfNeutronWallPointCollection, kTRUE);
  } else {
    FairRootManager::Instance()->RegisterAny("EtfNeutronWallPoint",
                                             fEtfNeutronWallPointCollection, kTRUE);
  }
  
  if ( ! gMC->IsMT() ) {
    FairRootManager::Instance()->Register("EtfNeutronWallFirstHits", "EtfNeutronWall",
                                          fEtfNeutronWallFirstHits, kTRUE);
  } else {
    FairRootManager::Instance()->RegisterAny("EtfNeutronWallFirstHits",
                                             fEtfNeutronWallFirstHits, kTRUE);
  }
  
}


TClonesArray* EtfNeutronWall::GetCollection(Int_t iColl) const
{
  if (iColl == 0) { return fEtfNeutronWallPointCollection; }
  else { return NULL; }
}

void EtfNeutronWall::Reset()
{
  fEtfNeutronWallPointCollection->Clear();
  fEtfNeutronWallFirstHits->Clear();
}

void EtfNeutronWall::ConstructGeometry()
{
	TString fileName=GetGeometryFileName();
	if (fileName.EndsWith(".geo")) {
		LOG(info)<<"Constructing geometry from ASCII file "<<fileName;
		ConstructASCIIGeometry();
	} else if (fileName.EndsWith(".root")) {
		LOG(info)<<"Constructing geometry from ROOT file "<<fileName;
		ConstructRootGeometry();
	} else {
		LOG(fatal) << "Geometry format not supported.";
	}

}

EtfNeutronWallPoint* EtfNeutronWall::AddHit(Int_t trackID, Int_t detID,Int_t id1, Int_t id2,
                            TVector3 posIn,TVector3 pos_out, TVector3 momIn,
                            TVector3 momOut,Double_t time, Double_t length,
                             Double_t eLoss)
{
  TClonesArray& clref = *fEtfNeutronWallPointCollection;
  Int_t size = clref.GetEntriesFast();
  return new(clref[size]) EtfNeutronWallPoint( trackID, detID, id1,  id2,
                             posIn, pos_out,  momIn,
                             momOut, time,  length,
                              eLoss);
}

EtfNeutronWallFirstHits* EtfNeutronWall::AddHit1(Double_t x0, Double_t y0, Double_t z0, Double_t t0,
                                   Double_t x1, Double_t y1, Double_t z1, Double_t t1,
                                   Double_t x2, Double_t y2, Double_t z2, Double_t t2,
                                   Double_t x3, Double_t y3, Double_t z3, Double_t t3,
                                   Double_t x4, Double_t y4, Double_t z4, Double_t t4,
                                   Double_t x5, Double_t y5, Double_t z5, Double_t t5)
{
  TClonesArray& clref = *fEtfNeutronWallFirstHits;
  Int_t size = clref.GetEntriesFast();
  return new(clref[size]) EtfNeutronWallFirstHits(x0,y0,z0,t0,x1,y1,z1,t1,x2,y2,z2,t2,x3,y3,z3,t3,x4,y4,z4,t4,x5,y5,z5,t5);
}

void EtfNeutronWall::Getpripoint()
{
  // <DB> dump of step info.
  static Int_t iStepN;

  // Particle being tracked
  const char *sParticle;
  switch(gMC->TrackPid()){
    case 2212:          sParticle="proton"    ;break;
    case 2112:          sParticle="neutron"   ;break;
    case 22:            sParticle="gamma"     ;break;
    case 50000050:      sParticle="ckov"      ;break;
    case 111:           sParticle="pi0"       ;break;  
    case 211:           sParticle="pi+"       ;break;  
    case -211:          sParticle="Pi-"       ;break;  
    case 1000010020:            sParticle="deuteron"        ;break;
    case 1000010030:            sParticle="triton"        ;break;
    case 1000020030:            sParticle="he3"        ;break;
    case 1000020040:            sParticle="alpha"        ;break;
    default:            sParticle="not known" ;break;
  }


  TString flag="-I- STEPINFO: tracking status: ";
  if(gMC->IsTrackAlive()) {
    if(gMC->IsTrackEntering())      flag="enters to";
    else if(gMC->IsTrackExiting())  flag="exits from";
    else if(gMC->IsTrackInside())   flag="inside";
  } else {
    if(gMC->IsTrackStop())          flag="stopped in";
  }
  
  //Where am i ?
  Int_t vid=0,copy=0;
  //current volume and his mother are always there
  TString path=gMC->CurrentVolName(); path.Prepend("-");path.Prepend(gMC->CurrentVolOffName(1));
  vid=gMC->CurrentVolOffID(2,copy);  if(vid) {path.Prepend("-");path.Prepend(gMC->VolName(vid));}
  vid=gMC->CurrentVolOffID(3,copy);  if(vid) {path.Prepend("-");path.Prepend(gMC->VolName(vid));}
  
//  printf("step  (nr=%i): %s (%i) %s %s m=%.6f GeV q=%.1f dEdX=%.4f Etot=%.4f \n",
//  iStepN,sParticle,gMC->TrackPid(),flag.Data(),path.Data(),gMC->TrackMass(),gMC->TrackCharge(),gMC->Edep()*1e9,gMC->Etot());
  
  Double_t gMcTrackPos[3]; gMC->TrackPosition(gMcTrackPos[0],gMcTrackPos[1],gMcTrackPos[2]);
  Double_t  gMcTrackPosLoc[3]; gMC->Gmtod(gMcTrackPos,gMcTrackPosLoc,1);
//  printf(" : track Position (MARS) x: %5.3lf, y: %5.3lf, z: %5.3lf (r: %5.3lf) ---> (LOC) x: %5.3f, y: %5.3f, z: %5.3f",
//          gMcTrackPos[0],gMcTrackPos[1],gMcTrackPos[2],
//          TMath::Sqrt(gMcTrackPos[0]*gMcTrackPos[0]+gMcTrackPos[1]*gMcTrackPos[1]+gMcTrackPos[2]*gMcTrackPos[2]),
//          gMcTrackPosLoc[0],gMcTrackPosLoc[1],gMcTrackPosLoc[2]);
  
//  printf("step (nr=%i): tid=%i flags alive=%i disap=%i enter=%i exit=%i inside=%i out=%i stop=%i new=%i",
//		 iStepN,  gMC->GetStack()->GetCurrentTrackNumber(),
//		 gMC->IsTrackAlive(), gMC->IsTrackDisappeared(),gMC->IsTrackEntering(), gMC->IsTrackExiting(),
//		 gMC->IsTrackInside(),gMC->IsTrackOut(),        gMC->IsTrackStop(),     gMC->IsNewTrack());
  
  Float_t a,z,den,rad,abs; a=z=den=rad=abs=-1;
  //Int_t mid=gMC->CurrentMaterial(a,z,den,rad,abs);
//  printf("step (nr=%i): mid=%i a=%7.2f z=%7.2f den=%9.4f rad=%9.2f abs=%9.2f\n\n",iStepN,mid,a,z,den,rad,abs);

/*  
  TArrayI proc;  gMC->StepProcesses(proc); 
  printf("-I- STEPINFO: Processes in this step:\n");
  for ( int i = 0 ; i < proc.GetSize(); i++)
  {
    printf("-I- Process involved --->   code: %i : name: %s\n", proc.At(i) , TMCProcessName[proc.At(i)]);
  }
  printf("-I- STEPINFO: end of process list\n");  
*/

  //Int_t Nprim=gMC->GetStack()->GetNprimary();

  Int_t trackNo=gMC->GetStack()->GetCurrentTrackNumber();
  Int_t motherNo=gMC->GetStack()->GetCurrentParentTrackNumber();
  Double_t time=gMC->TrackTime();

/*
  cout << "Prim: " << Nprim<<endl;
  cout << "TrackID: " << trackNo<<endl;
  cout << "MotherID: " << motherNo <<endl;
  cout << "time: " << time <<endl;
  cout << "PiD: "<<gMC->TrackPid()<<endl;
*/  

  if(gMC->TrackPid()==2112 && motherNo<0){
  // we have a primary neutron
     TArrayI proc;  gMC->StepProcesses(proc); 
     for ( int i = 0 ; i < proc.GetSize(); i++){
//        printf("-I- Process involved --->   code: %i : name: %s\n", proc.At(i) , TMCProcessName[proc.At(i)]);
     
        if(proc.At(i)!=22 && proc.At(i)!=23 && proc.At(i)!=31 && proc.At(i)!=43 &&  proc.At(i)!=13){
           LOG(INFO) << "new primary neutron interaction: " << proc.At(i) <<"  "<< TMCProcessName[proc.At(i)] << FairLogger::endl;
	
	}
         
// make histogram with first interaction
  
        if(proc.At(i)==22 && trackNo<6){
	   // elastic scattering
//           cout<<"primary neutron interaction elastic"<<endl;
           firstHitX[trackNo] = gMcTrackPos[0];
           firstHitY[trackNo] = gMcTrackPos[1];
           firstHitZ[trackNo] = gMcTrackPos[2];
           firstT[trackNo] = time;
  
        }
	
        if(proc.At(i)==13 && trackNo<6){
	   // elastic scattering
//           cout<<"primary neutron interaction hadronic"<<endl;
           firstHitX[trackNo] = gMcTrackPos[0];
           firstHitY[trackNo] = gMcTrackPos[1];
           firstHitZ[trackNo] = gMcTrackPos[2];
           firstT[trackNo] = time;
  
        }

        if(proc.At(i)==23 && trackNo<6){
	   // elastic scattering
//           cout<<"primary neutron interaction inelastic"<<endl;
           firstHitX[trackNo] = gMcTrackPos[0];
           firstHitY[trackNo] = gMcTrackPos[1];
           firstHitZ[trackNo] = gMcTrackPos[2];
           firstT[trackNo] = time;
  
        }
     }  
  }
  
  iStepN++;

}

void EtfNeutronWall::BeginEvent()
{
;
}
void EtfNeutronWall::FinishEvent()
{
;
}
void EtfNeutronWall::EndOfEvent()
{

  LOG(info) << "EtfNeutronWall: " << fEtfNeutronWallPointCollection->GetEntriesFast()
            << " points registered in this event";

  fEtfNeutronWallPointCollection->Clear();

}

void EtfNeutronWall::DefineSensitiveVolumes()
{
  TObjArray* volumes = gGeoManager->GetListOfVolumes();
  TIter next(volumes);
  TGeoVolume* volume;
  while ( ( volume = static_cast<TGeoVolume*>(next()) ) ) {
    if ( CheckIfSensitive(volume->GetName()) ) {
      LOG(debug2)<<"Sensitive Volume "<< volume->GetName();
      AddSensitiveVolume(volume);
    }
  }
}

Bool_t EtfNeutronWall::CheckIfSensitive(std::string name)
{
	cout<<"check if sensitive, for "<< name<<endl;
  if(TString(name).Contains("padle_h_box5")||TString(name).Contains("padle_h_box2"))
  {
		 cout<<"it is sensitive"<<endl;
    return kTRUE;
  }
	 cout<<"it is not sensitive"<<endl;
  return kFALSE;
}

ClassImp(EtfNeutronWall)
